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We study the deconfining phase transition in 3+1 dimensional pure SU{N) Yang-Mills theory 
. . . using a gauge invariant variational calculation. We generalize the variational ansatz of to mixed 

04 ' states (density matrices) and minimize the free energy. For > 3 we find a first order phase 

' transition with the transition temperature of Tc — 450Mev. Below Tc the Polyakov loop has 

vanishing expectation value, while above Tc, its average value is nonzero. According to the standard 
lore this corresponds to the deconfining transition. Within the accuracy of our approximation the 
j^j^' entropy of the system in the low temperature phase vanishes. The latent heat is not small but, 

I rather, is of the order of the nonperturbative vacuum energy. 



X 



I. INTRODUCTION 



(N 

• Attempts to understand the nature of the deconfining phase transition in QCD date back almost 25 years. Since the 

' pioneering work of Polyakov much effort has been made to study the basic physics as well as the quantitative 

■ characteristics of the transition. The high temperature phase of QCD is widely believed to resemble an almost free 

00 ' plasma of quarks and gluons. At asymptotically high temperatures this is confirmed by explicit perturbative calcu- 
lations of the free energy Perturbation theory in its simplest form, however, is valid only at unrealistically high 
temperatures. In recent years a different promising avenue has been explored. This incorporates analytical resum- 

, mation of the effects of the gluon screening mass into the 3D effective Lagrangian, which is then solved numerically 

f~| • by 3D lattice gauge theory methods & ^ . The results of this approach seem to be in good agreement with direct 

Qh! 4D lattice gauge theory calculations ^ all the way down to 2Tc. The numerical results indicate that, interestingly 

^ > enough, although asymptotically the free energy does approach that of the free partonic plasma, the deviations from 

I the Stefan-Boltzmann law even at temperatures of order lOTc are quite sizable, of order of 15%. One can inter- 

. pret this as an indication that the interesting physics of the transition region remains important even at these high 

L| ■ temperatures. 



Although we are quite advanced in the understanding of the high temperature phase, the transition region itself is 
' understood very poorly. This region of temperatures Tc < T < 2Tc, is of course the most interesting one, since 
5^ , it is in this region that the transition between "hadronic" and "partonic" degrees of freedom occurs. In fact, it is 
quite possible that the dynamics of the transition region is dominated by the magnetic vortices [Q , which are 
responsible for the abrupt increase in entropy . The study of the transition region is a complicated and inherently 
nonperturbative problem which has not been tackled so far by analytical methods. 

The purpose of this paper is to explore an analytical investigation of the phase transition region. Since the nature of 
the problem is nonperturbative, a nonperturbative method is called for. Unfortunately there are not many of those 
on the market. It does however seem to us that the variational method developed several years ago in H is quite 
suitable. We are encouraged by the fact that the variational approximation of manages to reproduce all known 
infrared features in lower dimensional confining theories |l^, 1^, ^ , and also gives a stable picture of instanton 
dynamics in QCD ||l6|| . Even though, as any variational method, it has an element of guesswork about it, it is a first 
principle calculation. For such a difficult physics problem as the one at hand, clearly even an exploratory investigation 
which provides a qualitative picture should be very valuable and is well worth attempting. 

In this paper, therefore, we will study the deconfining phase transition in a pure SU (N) Yang-Mills theory using the 
variational approach of |9j suitably extended to finite temperature. We will minimize the relevant thermodynamic 
potential at finite temperature, i.e. the Helmholtz free energy, on a set of gauge invariant density matrices. 



The results of our calculation are very interesting. We find that the phase transition in gluodynamics is closely related 
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to the symmetry breaking phase transition in the effective cr-model for the Polyakov loop variable U . This u-model 
arises explicitly in our calculation, and its parameters depend on the variational parameters of the trial density matrix. 
At low temperature the effective cr-model is in the disordered phase, and {U) = 0. The entropy of the Yang-Mills 
thermal density matrix vanishes in the low temperature phase. At the critical temperature of order ~ 450Mev, 
the best variational density matrix changes abruptly, and above Tc corresponds to the ordered phase of the cr-model. 
Thus at T > Tc we find {U) ^ 0. The transition itself turns out to be strongly first order for large N . This is 
consistent with recent lattice gauge theory calculations px| ]. These results are not surprising in themselves. They are 
expected from the general relation between the average value of the Polyakov loop and the (de) confining properties 
of the density matrix. Moreover, the order of the transition is predicted by the Svetitsky-Yaffe conjecture jl^. The 
vanishing of the entropy in the confining phase is also very natural, since the glueballs are heavy and their contribution 
to the entropy should be suppressed by the Boltzmann factor exp[— M^/T]. Nevertheless, we are not aware of any 
first principle analytic calculation in QCD which provides an explicit realization of these general arguments. 

This paper is structured as follows. In sec. 2 we recap the variational ansatz of for energy minimization at zero 
temperature and recall the methods of evaluating the expectation values in the variational wave function as well as 
the results of the minimization procedure. In sec. 3 we generalize the variational ansatz so that it incorporates not 
only pure states but also general density matrices. We introduce the two variational parameters with respect to which 
the free energy is to be minimized. In sec. 4 we recast the problem of the minimization of free energy in terms of an 
effective nonlinear cr-model. In sec. 5 we calculate the Helmholtz free energy of the trial density matrix and perform 
the minimization. We conclude in sec. 6 by discussing our results. 



II. THE VARIATIONAL CALCULATION AT ZERO TEMPERATURE 



The variational method proposed in is based on minimization of the expectation value of the Hamiltonian on a set 
of gauge invariant states. The set of states is chosen so that their functional form is a relatively simple generalization 
of the vacuum of the free theory. Ideally, one would like to choose a set of variational states that is both simple 
enough to perform explicit calculations and rich enough to allow for variations of its parameters to span interesting 
physics. One starts then with a Gaussian state 



{x){G-')lHx,v)A){y) \ , (1) 



where the set of functions G^^^^{x) are variational parameters. Expectation values of various operators are easily 
calculable in this set of states. One also expects that the freedom allowed by variation of G is wide enough to probe 
the nonperturbative physics of confinement. The problem is that these states are not gauge invariant and, therefore, 
as such do not belong to the physical Hilbert space of gluodynamics. To remedy this problem one projects the states 
onto the gauge invariant subspace by gauge transforming them and integrating over the whole gauge group: 



^'^{x)G-''^\x-y)A^\y)\ (2) 



with Af" defined as 



and 



AY'^{x)^S'^\x)A\{x) + \'i{x), (3) 



S'^^x) = \tr [t'^U^t^'U) ■ Xl{x) = -tr {t^^U^ d,U) (4) 

with T — traceless N x N hermitian matrices — normalized by 

tr(TV'') = 25"^ (5) 
and obeying the completeness condition for SU{N) 

1 



2{Sii6jk - j^^ijSki) ■ (6) 
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The integration in eq. (H) is performed over the space of special unitary matrices with the SU {N) group invariant 
measure. This integration projects the original Gaussian state onto a colour singlet. Due to the projection operation, 
the calculation of expectation values in this state is much more involved than in the case of a simple gaussian. Thus 
it is a very difficult task to perform the full fimctional minimization with respect to Gj'j,(a:;). This can be done and 
was indeed performed in simpler calculations in 2+1 dimensional models Jl2| , [l^ . In these models the approach was 
also extended to provide for the calculation of the confining potential between the external charges [Q. It has also 
been applied to 1+1 dimensional models with fermions jl5|. In these simple models the variational approximation 
has performed very well, reproducing all known results. 

In nonabelian theories, unfortunately, full functional minimization is beyond our calculational abilities. Therefore to 
make headway, the width of the Gaussian G was restricted in |^ to a one parameter form as follows 

Gt^ix~y)^S-'S,,Gix-y) (7) 

with the Fourier transform of G taken to be 

^ ^^>-\M \ie<M^ ■ 

The dimensional parameter M is thus the only variational parameter in the calculation. 

The expectation value of any gauge invariant operator O in the variational state eq. (H) is given by the following 
integral 

(O) = ^ [ DU{0)a , (9) 



z ^ 

^ j UAe^^ / <ixrfyA^"(a;)G-i(a;-y)Af"(y)Qg-i / dx' dy' A)(x')G-^ (x' -y')A)(y') _ 

Note that only one group integral is present. This is because for a gauge invariant operator O only one of the states 
has to be gauge projected. 

Since the gauge transform of a vector potential is a linear function of A eq. , for fixed U {x) this is a Gaussian 
integral, and can therefore be performed explicitly for any reasonable operator O. The nontrivial part of the calculation 
is the path integral over the group variable U{x). Consider first the normalization factor Z. After integrating over 
the vector potential we obtain: 

Z = J DU cxp{-S[U]} (10) 



with the action 



Here 



S[U]^^Tr\nM + ^X[G + SGS'^]-^X. (11) 



St^'ix, y) = S'^\x)5,jS{x - y), M^'^x, y) = [S' '''{x)S^'{y) + 5'^']G-'{x - y)5,, , (12) 

and we use matrix product notation, which implies summation over the colour and Lorentz indices, as well as integra- 
tion over the spatial coordinates. In the rest of this paper Tr stands for summation over all indices, including the space 
coordinates, while tr stands for summation over the colour indices alone. The path integral eq. (|lC| ) defines a partition 
function of a nonlinear a- model with the target space SU {N)/Zj^ in three dimensional Euclidean space. The fact that 
the target space is SU{N)/Zn rather than SU{N), follows from the observation that the action eq. ( |ll| ) is invariant 
under local transformations belonging to the centre of SU (N) . This can be trivially traced back to invariance of Af 
under gauge transformations that belong to the centre of the gauge group. 

The matrix U {x) has a well defined gauge invariant meaning. It arises in eq. (|lO|) as a relative gauge transformation 
between the bra and the ket gaussian wave functionals. As such it plays exactly the same role as Polyakov's loop P 
at finite temperature. We will therefore refer to it as the Polyakov loop in the following |l^. The contribution of a 
given U {x) to the partition function eq. (|lO|) and to other expectation values corresponds to the contribution due to 
the off diagonal matrix element between the initial Gaussian and the Gaussian gauge rotated by U{x). Therefore, if 
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matrices U{x) which are far from unity give a significant contribution to the partition function, it means that the off 
diagonal contributions are large and gauge projection is indeed mandatory to get physically sensible results. 

The action of this cr-model is rather complicated. It is a nonlocal and a nonpolynomial functional of U (x) . The 
following set of approximations was employed in ^ to estimate this path integral. First, since the bare coupling 
constant of the Yang Mills theory is small and it also enters as a coupling constant of the nonlinear a-model, the high 
momentum modes of the field U are integrated out perturbatively. The result of the perturbative one loop integration 
of the modes with spatial momenta fc > M is a much simpler u-model which involves low momentum modes only. 
The action of this cr-model is 



M 



2g^M) 



tr / d^x d^U^x)diU{x) 



(13) 



where M also serves as the ultraviolet cutoff. To be more precise, due to the Zn local symmetry of the original theory 
eq. ([Ill), the action for the low momentum modes is slightly different. The derivatives should be understood as Zn 
covariant derivatives. The most convenient way to write this action, would be to understand U{x) as belonging to 
U{N) rather than SU{N), and to introduce a C/(l) gauge field by 



1 M 
Sl = 7;—tr 

2 ff^ 



d^x {d, ~ iAi)U\x){d, + iA,)U{x) . 



(14) 



This defines a cr-model on the target space U{N)/U{1), which is isomorphic to SU{N)/Zj^. However, since the 
difference between the two actions is not important at large number of colours, we will disregard this subtlety in the 
following and use eq. 



As a result of the running of the QCD coupling constant, the effective coupling 5^(M) is given by the one loop QCD 
formula [p^ . Clearly for the perturbative integration to be consistent we should find that the energy is eventually 
minimized for high enough value of the scale M, so that as{M) ^ 1. The actual value found in ^ is as{M) ~ .25, 
which is indeed reasonably small. To solve the resulting low energy cr-model the mean field approximation and 
perturbation theory are used. We will describe this procedure here in some detail, since we will use it in later sections 
to analyze the finite temperature case. 

The Hamiltonian of the pure Yang-Mills theory is 



d-^x 



-e: 



a2 



-B 



a2 



(15) 



where 



5A-{x) 



-e,,k{d,Al{x) - dkA'^^{x)+gr'''A){x)Al{x)} . 



(16) 



To calculate the expectation value of H in the variational state we first perform the integral over the gauge fields A 
at fixed value of the Polyakov loop U . The result of this integration (for details see ^) for the energy density is 

(2H) _ 3(iV^^_,^^^^^^^^,_^^^,^,^^^_^^l^^^ 



V 2 

- ±. I d''xd^y{\'iH{x)G-\x - y)XlH{y))u + \ 



d^x{\1^(x)\Ux))u , (17) 



where the averaging over the U field should be performed with the cr-model action, eq. (|l^). As the high momentum 
modes are only considered to one loop order, the cr-model action reduces to 

S^\j dxdy\lH{x)G-\x - y)KH{y) + ^^5^^^ / '^'^ d,ul{x)d.UL{x) . (18) 
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The subscripts H and L indicate that the corresponding expressions contain only high and low momentum components 
of the fields respectively. Integration over the high momentum modes gives 



2(H) 
V 



N 



1 kt"^ 



(19) 



where now only the contribution of the low momentum modes remains to be evaluated. 



This last step of the calculation is the most difficult one. The nonlinear cr-model eq. (13) possesses a global SUl{N) (g) 
SUr{N) symmetry. Thus it has a symmetry restoring phase transition as a function of M. The coupling constant 
g^iM) enters this action in the same way as temperature enters in the classical statistical problem of the SU{N) 
spin system. High values of M correspond to low values of g'^{M) and thus to the low temperature regime in the 
cr-model. In this regime the global SUL{N)(S)SUuiN) symmetry is spontaneously broken down to the vector subgroup 
SUv{N). The model itself is weakly interacting and thus one can use perturbation theory. However, it is well known 
that the model undergoes a phase transition from the ordered low temperature phase to the disordered phase at some 
critical temperature g^. The disordered phase is clearly nonperturbative. To treat it quantitatively the mean field 
approximation was used in . The partition function of the cr-model is rewritten by introducing a (hermitian matrix) 
auxiliary field a which imposes a unitarity constraint on U{x) 



(20) 



S[U,a] = 



M 



2g^{M) 



—tr / d?'x[d^U\x)d^U{x)+a{U^U ~\)] . 



In this step we have disregarded the fact that U is an SU{N) matrix rather than U{N). Thus, strictly speaking, our 
calculation applies only at large N. The mean field equations for this action are 



(C/tC/) 



(21) 



(aU) = 0. 



(22) 



From eq. ( p^ ) it follows that either (a) = 0, (U) ^ (the ordered, broken symmetry phase with massless Goldstone 
bosons), or (a) ^ 0, {U) — (the disordered, unbroken phase with massive excitations). 



In the disordered phase the expectation value of a is proportional to the unit matrix 
Eq. ( ^ ) then becomes 



27V- 



:5'(M) 



M 



1 



N'^g'^(M) 



M Jo (27r)3 + 



a M 

1 arctan — 

M a 



N. 



(23) 



(24) 



The gap equation, eq. (p^), has a solution only for couplings g'^(M) larger than the critical coupling g"^, which is 
determined by the condition that a = 



' 47r 4 iV 

The low momentum mode contribution to the ground state energy is 



(2^)3 A:2 + ^: 



2^2 



rM 



1 



M' 



-M'^ - a^M + a-^ arctan — 



(25) 



(26) 



The final mean field expression for the ground state energy density is 



4^ 



15 NP as{M) 



(27) 



6 



where as{M) is the QCD couphng at the scale M, is given by eq. (p5[), and a is determined by 

— arctan — = ^ '^^ " . 28) 

This function has the minimum at the critical point q;s(M) — a^. With the one loop Yang-Mills function and 
Aqcd = 150Mev, we find for iV ^ 3 

24 

M, = Agczje- ^ 8.86Aqci3 = 1.33Gev . (29) 

Eq. (|2^) is valid only in the disordered phase. Here we expect the mean field approximation to be reasonably good, 
since the spectrum of the cr-model is massive and fluctuations arc not expected to be overwhelmingly large. 

In the ordered phase, on the other hand, the situation is different. The spectrum contains 0{N^) "Goldstone bosons" . 
In the mean field approximation in addition it also has 0{N'^) massive particles — the analogs of S particles in QCD. 
In the mean field the 5's are all stable. In reality of course, since there is nothing to prevent them from decaying into 
the Goldstone bosons, they are not at all stable. In fact, the further into the ordered phase one gets, the less stable 
they are due to the increase of the phase space available for decay. Thus we expect the mean field to be unreliable in 
the ordered phase. Asymptotically at very large M the mass of 6 becomes much higher than the ultraviolet cutoff of 
the cr-model. For these values of M the massive modes decouple even within the mean field approximation, and one 
recovers the perturbative result. This decoupling, however, occurs extremely slowly. We have checked that the mean 
field formulae tend to the perturbative ones only when the mass of the would be stable massive states is about 1000 
times the ultraviolet cutoff. This is, of course, unphysical, confirming our suspicion that the mean field approximation 
cannot be used on the ordered side of the transition. 

The simplest option to calculate the contribution of the low momentum modes to the energy in the ordered phase is 
to use perturbation theory. The raison d'etre for this is the following. Perturbation theory is certainly appropriate 
for large enough values of M, where the expectation value of the U field is of order unity. On the other hand, 
we know from the numerical studies of [p^ that the phase transition in our cr-model is strongly first order. More 
specifically, according to ||T^ the transition occurs when the expectation value of U is greater than .5. We thus expect 
perturbation theory to be qualitatively reliable all the way down to the transition point. Since we have neglected 
previously many perturbative contributions of order g^, in order to be consistent we have to limit ourselves to the 
leading order contribution only. Calculating with this precision the energy of the trial states for M > Mc we get 



E(M) _ -I 



. (30) 
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Thus for M > Mc the energy is a monotonically increasing function of M. 

Eq. ( p^ ) is thus the result of the minimization of the energy over the whole range of the variational parameter M . 

III. THE VARIATIONAL ANSATZ FOR THE DENSITY MATRIX 

In order to extend the variational analysis to finite temperature we have to generalize our ansatz so that it includes 
mixed states. In scalar theories the Gaussian approximation has a long history of applications at finite temperature 
[pO| , We generalize our ansatz along the same lines. 

We start by considering the density matrices which in the field basis have Gaussian matrix elements 

g[A,A'] = expj - ^ ^A^ix)G-^'^\x,y)A^iy) + A[^x)G-^'^'{x,y)A'^\y) - 2A^ix)H^^ix,y)A'^'{y)^ ■ (31) 

As before, we take the variational functions diagonal in both colour and Lorentz indices, and translationally invariant 

G-^""'{x,y) = S''%,G-\x - y) , (32) 
Hf^{x,y) = S^%,H{x-y). 
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Then 

~q[A, A'] = exp I - i j {AG-^A + A'G'^A' - 2AHA') | . (33) 

Note that for H = this density matrix represents a pure state, since it can be written in the form 

g^\'^[A]><^[A]\ (34) 

with "^[A] a gaussian wave function. At nonzero H the density matrix is, however, mixed. The magnitude of H, 
therefore, determines the entropy of this trial density matrix. 

We now make an additional simplification in our ansatz. First, we restrict the functions G~^{x) to the same functional 
form as at zero temperature, eq. (^). Further, we will take H{k) to be small and nonvanishing only at low momenta 

^^''>-\H<^M ifP<M2 ■ ^^^^ 



The logic behind this choice of ansatz is the following. At finite temperature we expect H{k) to be roughly proportional 
to the Bolzmann factor exp{— E(/c)/3}. In our ansatz, the role of one particle energy is played by the variational function 
G^^{k). We will be interested only in temperatures close to the phase transition, and those we anticipate to be small, 
Tc < M. For those temperatures one particle modes with momenta k > M are not populated, and we thus put 
H(k) = 0. For k < M the Bolzmann factor is nonvanishing, but small. Further, it depends only very weakly on the 
value of the momentum. We will have, of course, to verify a posteriori that our assumptions about the smallness of 
Tc and H are justified. As we will see later, this turns out indeed to be the case with reasonable accuracy. 



As before, we explicitly impose gauge invariance by projecting g onto a gauge invariant sector 

g[A,A'] = J PC/Wexp| - A^' G^^A^' + A'^" G^^A'^" - 2A^'iJA'^'"| 



(36) 



One of the group integrations in eq. ( pq ) is redundant, since we will only calculate the quantities of the form TrgO, 
with O being gauge invariant. We thus have the following ansatz for the density matrix 



g[A,A'] 



VUexp 



AG-^A + A'^G-^A'" - 2AHA'^ 



(37) 



This expression is not explicitly normalized to unity. Nevertheless, we find it convenient to refer to it as density 
matrix while explicitly inserting a normalization factor whenever necessary. Thus the average of a gauge invariant 
operator O is given by 



(O) 



A,U 



= Z-^ J VUVA 0{A, A') exp I - i ^ AG'^A + A'^G^M'^ - 2AHA'^^ 

where Z is the normalization of the trial density matrix g, i.e. 

Z = ^^Q = j VUVAexp | ^ ^ ^ AG'^A + A^G~^A^ - 2AHA^^ 



VUVA exp 



AAA + \{G-^ - Lo/\-'^uj^)\ 



with 



A = A + XluA~ 



HG, 



A = 2G-^[1-—{S + S 
CO = {G-^S~H). 



A'=A 



(38) 



(39) 

(40) 
(41) 
(42) 
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The A integration can then be easily performed to yield 

„ ?! 

(43) 



Tig^ J PC/expj - iA(G-i -wA-iw'^)A- ^Trln^l 



We now adopt the same strategy for treating the high momentum modes of J7 as at T = 0. Namely, they are integrated 
perturbatively to one loop accuracy. The result is the effective a-model for the matrices U with momenta below M. 
The coupling constant g of this cr-model gets renormalized as before according to the one loop Yang- Mills /3-function, 
and thus has to be understood as g{M). Additionally, due to independence of H on momentum, for low momentum 
modes of U the function H{x — y) is equivalent to H6^{x — y). 

The final approximation has to do with the fact that H is assumed to be small. For arbitrarily large H the variational 
calculation is forbiddingly complicated even with all the above mentioned simplifications. This is because the gauge 
projection renders the calculation of entropy in the general case unfeasible. However, at small H we only need to 
calculate the leading term in entropy. This calculation can indeed be done, and is described in the next section. Since 
we are only calculating the leading order contribution in H, we only have to consider corrections to the cr-model 
action of first order in H. With this in mind, the normalization factor becomes 

TTg= fvUexp(-^\(^^ + ^{S + S^)y+^HGtiiS + S^)\. (44) 



IV. THE EFFECTIVE a-MODEL 

The normalization Z can be interpreted as the generating functional for a theory defined by the action S{U) 

Z = Trg = J VUe-^^^^ , (45) 



where 



Af 1 

Saj) = — AA + -A H{S + S'^)X -HM^tiS . (46) 

4 8 47r2 



We simplify this expression using 



AA ^tr(9C/aC/t) , (47) 
9 

XS^X = XSX = --^tr\(U''dU - dU'^UMdUU'^ - UdU'^)] , (48) 
2.9 L -I 

trS* ^ trS*^ = tiUhiU - 1 . (49) 



Inserting these into the action we get 

S{U) = ^tridUdU'') - ^ti-\(U^dU - dU''U)(dUUl - UdW^)] - -K^HM^trUhrU , (50) 
2(7^ 8g^ L J 47r^ 

where U independent pieces have been dropped. 

As noted before, this effective a-model is nothing but the effective theory for the low momentum modes of the Polyakov 
loop variable. Its status and applicability region are different from the usual perturbative effective actions, see e.g. 
|P2| . The standard effective action is calculated in perturbation theory and is valid at high temperature. Our effective 
action eq. (|5^) depends on the variational parameters M and H , and in a sense is a variational effective action. Also 
due to our restrictions to small values of H, a priori we do not expect it to be valid at high temperatures but, it 
rather, should represent correctly the physics in the phase transition region. 

Another important difference is that our effective tr-model does not have the local gauge invariance U{x) — s- 
V'' {x)U {x)V (x) which is usually associated with the effective action for the Polyakov loop. The reason for this 
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is that our setup is different from tfiat of the standard finite temperature calculation. The way this gauge invariance 
usually appears is the following. Consider the calculation of any gauge invariant observable 

{O) - I VUTr[eM-mOg{U)] , (51) 

where g{U) is the quantum mechanical operator of the gauge transformation represented by the matrix U. This 
expression for fixed U can be compared to the same expression but with U gauge transformed 

Tv[eM-mOg{V^UV)] = Tr[exp{-(3H}Og{V^)g{U)g{V)] = Tr[exp{-/3H}05({7)] . (52) 

The last equality here follows from the fact that both O and exp{/3H} are gauge invariant, and thus the operator 
(7(y^) can be commuted all the way to the left. The only effect of the transformation is then to change the basis over 
which the trace is being taken, which obviously leaves the trace invariant. 

Our setup is somewhat different. Expectation values are calculated as 

J VUTT[ggiU)0] (53) 

with g defined in eq. (^3|). This expression is altogether gauge invariant, since the integral over U correctly projects 
only the contribution of gauge singlet states. However the operator g is not itself explicitly gauge invariant. For that 
reason the gauge transformation operator giV'^) cannot be commuted through it, and thus 

Tt[gOg{V^UV)] ^ Tr[gOgiU)] (54) 

even for gauge invariant operators O. This manifests itself as absence of local gauge invariance in the action of the 
effective c-model, eq. (^. 

Nevertheless, we stress again that the meaning of the SU{N) valued field U is precisely the same as that of the 
Polyakov loop. 



V. THE CALCULATION OF THE FREE ENERGY 



To find the best variational density matrix we have to minimize the free energy with respect to the variational 
parameters M and H. The Helmholtz free energy F of the density matrix g is given by 

F=(H)-TS, (55) 

where H is the standard Yang-Mills Hamiltonian eq. ([Tsl), S is the entropy, and T is the temperature. 

Thus 

F = ^(TT{E^g) + Tr{B^g)j + T-Tr{g\ng). (56) 

First of all we need to perform the integration over the gauge fields, and reduce this expression to the average of a 
[/-dependent operator in the effective tr-model. In fact, as we shall see soon, to leading order in H the only nontrivial 
calculation we need to perform is that of the entropy. 

We will calculate the entropy up to the first nontrivial order in H . As we now show, the leading term at small H is 
0{HlnH). 

Let us denote by go the density matrix of the pure state with H = 0: 

eo = |0)(0|. (57) 

Here |0) does not denote necessarily the actual ground state, but rather a projected Gaussian state with arbitrary M. 
Now, since the matrix elements of the density matrix can be expanded in powers oi H , to leading order we can write 



g^ go + 5g, 



(58) 
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where Sg is 0{H). 

Imagine that we have diagonalized g. It will have one large eigenvalue ao = 1 — 0{H), which corresponds to the 
eigenstate 

|0') = |0)+O(H). (59) 
All the rest of the eigenvalues are at most 0{H). Then the entropy can be written as 

S = — Tr(^ In = — ao In ao — c^i In on . (60) 

i=l 

The second term is 0{H In H), and it is the cocfRcicnt of this term that wc will now calculate. Neglecting 0{H) 
corrections, we can substitute aj = H/M under the logarithm. Thus to leading logarithmic order 

S = -^ailnif/M. (61) 

i 

So that all we have to calculate is a^. Let 

\0') = \0)+H\x). (62) 

Then 

OO 

e = ao|0')(0'| + ^a,|a;,)(a;,| (63) 

i=l 

with 

(a;i|0')=0. (64) 

Note that {0\x) ^ 0, but 

(0|a;) + (a;|0) = , (65) 

since |0') has to be normalized at 0{H). Also 

{xi\0) +H{xi\x) = 0. (66) 

Thus the overlap (a:;j|0) is 0{H), and we have 

g = ao|0)(0| + h(\0){x\ + |a;)(0|) + ai\xi){xi\ . (67) 

Multiplying this by go we get 

go-g = ao-go + H\0){x\ + i?(0|a;)|0)(0| , (68) 
g-go = ao- go + H\x}{0\ + H{x\0)\0){0\. 

Thus, 

go - g + g- go - g = ao - go- ai\xi){xi\ . (69) 
Multiplying by go again, we get rid of |a;,)(a;i| to 0{H) 

ao - go = go - g + go - Q - Qo- Qo- Q = Qo- Q - Qo- (70) 

Then, 

ao = Tr(eo • g) ■ (71) 
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Since Trg = 1 we have 



which, inserted into eq. (|60|), gives 



^ Qfj = 1 - ao Tr(£io(l - O)) 



S = -(1 - Tr(go ■ Q))\nH/M . 



(72) 



(73) 



The derivation has been given for the normahzed density matrices Qo and q. In terms of our Gaussian matrices we 
should restore the normahzation factors Z and Zo, so that finally we have 

• Tr(£»o • q) 



It is easy to check that to 0{H) 



Then 



From eq. (H), it is clear that 



Tr£)o ■ Trgi 
Tr(£»o • q) = (Trgo)^ 



1 In i7/M. 



S = 



(Ir 



Trp 



Trp = 



'(4772 



l + fl"(^M^tr5- -AS'A 



Trpo 



Using eqs. (^ 49) we finally get 
S = 



< -^tTOJ^dU - dU^U)(dUU^ - UdU'') + -l-M^UrUhrU - 1) >u 
8.g^ 47r^ 



HItlH/M. 



(74) 

(75) 
(76) 

(77) 
(78) 



To leading order in H, the averaging over U in this expression has to be performed with the cr-model action at 
vanishing H. 

Expression eq. (|78| ) has the following striking property. For AI < it vanishes identically. The reason is very 
simple. The first term in eq. ( [78| ) is the product of the left handed SU{N) current and the right handed SU{N) 
current in the cr-model. Thus it transforms as an adjoint representation under each one of the SU{N) factors of the 
SULiN) ig) SUji{N) transformation. The same is also true for the second term in eq. ([78|). The cr-model action itself 
at 7? = obviously is invariant under the whole SUl{N) (E)SUr{N) group. Now, at M < Mc, the symmetry group is 
not spontaneously broken, and thus any operator which is not a scalar has a vanishing expectation value. It follows 
immediately that the entropy has an 0{H\nH/M) contribution only for M > Mc, when the SUl{N) (g) SUji{N) 
group is spontaneously broken down to SUviN). 

This observation makes our task considerably simpler. Since for M < Mc the entropy is zero, we do not have to 
consider at all the disordered phase of the effective cr-model. In this disordered phase the free energy coincides with 
energy, and thus the calculation is identical to the calculation at zero temperature. There is one slight subtlety here. 
We expect of course that since the energy alone must have a minimum on a pure state, the nonvanishing H should 
give always a positive contribution to the energy. Thus in the absence of the entropy contribution, H should vanish 
in the minimal energy state. Nevertheless since in the energy at T = was only minimized with respect to AI, we 
have to check this property explicitly. We have indeed performed this check and found that for small H the derivative 
of the expectation value of the energy with respect to H is strictly positive. 

Thus we only need to consider the effective cr-model in the ordered phase. As at T = we perform the calculations 
in the ordered phase to leading order in as- Since there are no 0{H\nH/M) corrections to energy at this order, 
the result for the energy in the disordered phase is identical to the result at zero temperature, eq. (^o|). Thus our 
expression for the free energy is 

Ar2 _ 1 / 1 1 \ 

F = — — — Af4 + T (— tr([/^ac/ - dU^UMdUU^ - UdU"^) + —AI^(ti-UhiU - l))u\H\nH/M . (79) 



We now average over U in the leading order perturbation theory. 
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VI. THE cr-MODEL PERTURBATION THEORY — MINIMIZATION OF THE FREE ENERGY AND 

THE DEBYE MASS 

We parameterize the U matrices as 



U 



e^p[-'-rr^}. (80) 



Although we only need the leading order, it is instructive to check that the order term in the expansion is indeed 
small. To this order we have 

U~(l- -(/)V° - i^^/r'^r'' + A^"0''^VV''t'=) . (81) 
V 2 8 48 / 

So that the cr-model action becomes 



^ = TT^tridUdU^) (82) 



M 
2^ 

4^^^^^+ 1929 



T T T T ~ T T T T 



The propagator of the phase field (j> is thus 



S"" . (83) 



To get the idea of the quality of this perturbative expansion we can calculate for example (S) . In this calculation one 
has to take into account the fact that the measure in the path integral over the phase is not the simple but 
rather the group invariant U (N) measure fi . To first order in it is 

^, = Vrexp{^ I d'xq^'ix)}. (84) 

Taking this into account we find that (S) gets no correction of order g^. We thus feel confident that the use of the 
perturbation theory in the ordered phase of the cr-model is an admissible approximation. In the following we will only 
keep leading order expressions. 

Calculating to leading order the entropy eq. ( [78| ) and keeping only the 0{N'^) terms we find 

(S)=-£M^i/ln§. (85) 



Introducing the dimensionless quantity 



we can write the expression for the free energy as 



(86) 



F = (H)-T(S) 

We now have to minimize this expression with respect to h and M . It is convenient to first perform the minimization 
with respect to h at fixed M. This obviously gives 
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Thus as a function of M only, the free energy becomes 



1207r2' ' e 67r 

Now minimizing with respect to M we find 



I^^O^M^i^. (90) 
dM e ^ ^ 



Thus for M > Mc the free energy of the best variational density matrix as a function of temperature is 

(91) 



'M>Ma — ^ ' ' 



3607r2 \ e J 

We now have to compare this value with the free energy for M < Mc- As we have discussed above, this is given by 
the expectation value of the Hamiltonian alone, and is minimized at M = Mc- Its value is 

---KKM.^-^Mt. (92) 



Comparing the two expressions we find 



Using the value of from eq. (29) we have 



307r2' 



1246 

Tc = — . (93) 
15 



Tc = 450Mev . (94) 



For T < Tc the free energy is minimized in the variational state with M = Mc- In our approximation this state is the 
same as at zero temperature. Its entropy vanishes, and the effective cj-modcl is in the disordered phase. The Polyakov 
loop vanishes, ([/) = and according to the standard wisdom this is a confining state. 



For T > Tc the best variational state is very different. The entropy of this state is nonzero, 

(95) 



iV2 /15T\3 



GTT-^e V e / 

The Polyakov loop is nonzero {U) ^ and thus the high temperature density matrix describes a deconfined phase. 

Finally, we note that in the deconfined phase our best variational density matrix has a nonvanishing "electric screening" 
or "Debye" mass. The Debye mass is conveniently defined as the "mass" of the phase of the Polyakov loop. This 
mass is nonvanishing in our calculation for the following reason. As long as iJ = 0, the effective cr-model action has a 
global SUl(N) (g) SUr{N) symmetry. Thus in the ordered phase of the a model the phases 0° are massless. However, 



as discusses above, the terms of order H in eq.(|50[) break this symmetry explicitly down to the diagonal SUviN). As 
a result the would be "Goldstone" phases 0" acquire mass. To calculate this mass it is convenient first to note that 
to Oig^) 

triU^dU ~ dU^U){dUUl ~ UdU'^) = -AirdU'^dU - i<?!>"</."9/5(/)''tr(r''TV"r'^ - t^t'tV'') . (96) 



The contribution of the SUl{N) (g) SUii{N) term to the mass cancels against the contribution of the measure eq.(p4[). 
Using eqs.(|50|, |9^ ) we then find to 0{g'^) and to leading order in H 

Ml^^a,{M)NMH . (97) 



As a function of temperature we have 



Ml^as(-T)N^-^T\ (98) 



14 



VII. DISCUSSION 



Let us summarize the results of our variational calculation. We find the phase transition at a temperature of about 
Tc ~ 450Mev. The transition is strongly first order at large N. The latent heat is AE = ■^:^{^^) ■ Below the 
transition the entropy is zero, the best variational state is the same as at zero temperature, and the average value of 
the Polyakov loop is zero. Above the transition, the entropy is nonzero and proportional to the number of "coloured" 
degrees of freedom, S oc N'^. The average value of the Polyakov loop is nonzero and the phase is deconfined. 

It is quite interesting that at high temperature our formulae numerically are quite close to the predictions of free 
gluon plasma. In particular, our value for the free energy, eq. (^), should be compared to the free gluon plasma 
expression 



The ratio between the two is 



The ratio of the entropies is the same. 



d^~ 0.85. (100) 



Interestingly we get the same ratio comparing our value for the Debye mass eq. (^8|) with the leading order perturbative 

one, M^^rt = fa,NT\ 

A,f2 

(101) 

The pressure approaches its asymptotic value according to the simple formula 



asympt 



One has to take the comparison eqs.(100 101 ) with a grain of salt. As explained above, our calculations were performed 
assuming small H. A priori we expect that this restriction should confine us to not too large temperatures. Interest- 
ingly, however, the minimization of the free energy resulted in the value H/M — 1/e independently of temperature. 



Thus, we feel ourselves justified to consider the comparison eq.(lOG) as meaningful 



The main features of these results are indeed what we expect from the deconfinement phase transition on general 
grounds. It is extremely gratifying that a simple minded calculation such as ours does qualitatively so well in such 
a complicated problem. It therefore appears that the projection of the trial density matrix on the gauge invariant 
Hilbert space is, just like at zero temperature, the crucial feature that dictates most if not all the important aspects 
of the low energy and low temperature physics. In the context of the present calculation the most important effect of 
the gauge projection is obviously vanishing of the entropy in the low temperature phase. We stress that this feature 
was not at all built into our initial ansatz, but followed naturally and unavoidably in the disordered phase of the 
effective u-model. 

Quantitatively, our calculation of course should be taken for what it is — an approximate implementation of the 
variational principle. The projection over the gauge group, which as we saw is so physically important, is what makes 
the calculational task difficult. The most severe simplifications that we had to impose are the perturbation theory in 
the ordered phase of the cr-model and the assumption of smallness of H. 

We believe that both these approximations affected the quality of our results. In particular, in the leading order of 
the perturbation theory the expectation value of the Polyakov loop U is equal to unity. The actual value of U on 
the ordered side of the transition according to [|l9j is close to one half. Thus our perturbative calculation is rather 
more reliable somewhat further away from the transition. In line with this we expect that the estimate for the critical 
temperature we obtained here is somewhat higher than we would get, had we treated the cr-model more accurately in 
the transition region. This is consistent with the fact that our result for Tc is by about 50% higher than the lattice 
value of 270Mev. 
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The smallness of h is also quite important. The value oi h = 1/e that we obtain is in fact a reasonably small number, 
so omitting the corrections in powers of h is fairly safe. On the other hand, the terms linear in h but not enhanced 
by In h, which we have ignored in the present calculation, have to be accounted for in a more careful way. Once these 
terms are taken care of our calculation can be extended to high temperatures. 

We believe that with some effort both these limitations can be overcome, at least to some extent. It should be possible 
to treat the cr- model action in a better way, perhaps along the lines of a continuum version of [p^ . As for terms 
linear in h, those should be accessible in perturbative expansion in h, once the nonanalytic 0(/iln/i) terms have been 
understood. We plan to continue investigations along these lines and hope to improve on the quality of our results. 

Notwithstanding this critique, we are very much encouraged by the results of this exploratory investigation. We believe 
that a reasonably simple improvement on the calculational methods has a good chance to significantly improve the 
results and bring them into a better quantitative agreement with the lattice data. 
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